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The ability to store information is of fundamental importance to any computer, be it classical or 
quantum. Identifying systems for quantum memories which rely, analogously to classical memories, 
on passive error protection ('self-correction') is of greatest interest in quantum information science. 
While systems with topological ground states have been considered to be promising candidates, a 
large class of them was recently proven unstable against thermal fluctuations. Here, we propose 
new two-dimensional (2D) spin models unaffected by this result. Specifically, we introduce repulsive 
long-range interactions in the toric code and establish a memory lifetime polynomially increasing 
with the system size. This remarkable stability is shown to originate directly from the repulsive 
long-range nature of the interactions. We study the time dynamics of the quantum memory in terms 
of difi'using anyons and support all our analytical results with extensive numerical simulations. Our 
findings demonstrate that self-correcting quantum memories can exist in 2D at finite temperatures. 



I. INTRODUCTION 

Quantum computers cannot be realized without the 
help of error correctioiP'. By encoding quantum informa- 
tion into logical states and designing correction circuits 
working on them, computations and information can in 
principle be protected from decoherence. However, the 
need of such an active control mechanism poses a major 
challenge for any physical implementation. It is therefore 
of greatest interest to look for passively protected systems 
which are intrinsically stable against the destructive in- 
fluence of a thermal environment. For this reason, the 
idea to encode quantum information in a topologically 
ordered ground state 1 5'o) of a suitable Hamiltonian has 
attracted a lot of intereslPf^. 

Important candidates a mon g such topological models 
are stabilizer HamiltonianJ^^, which are given by a sum 
of mutually commuting many-body Pauli operators. The 
advantage of such Hamiltonians is that the full energy 
spectrum is known and error correction schemes are read- 
ily derivecpli^. However, very recent results® ^ show that 
in one and two spatial dimensions no stabilizer Hamil- 
tonian with finite-ra nge in teractions (including the well- 
known Kitaev modeP^^EEl) can serve as a self-correcting 
quantum memory due to the errors induced by a thermal 
environment. In other words, increasing the size of such 
a system does not prolong the protection of its ground- 
state space from decoherence. These negative results 
point towards the fundamental question whether topolog- 
ically ordered quantum states, and hence self-correcting 
quantum memories, can exist at all on a macroscopic 
scale. In the following, we will demonstrate that self- 
correcting properties of 2D stabilizer Hamiltonians can 
indeed be established when we allow for long-range re- 
pulsive interactions between the elementary excitations 
(anyons). While the purpose of the present work is of 
principal nature, we note that such interacting models 
can be expected to be realized in physical systems. We 
discuss this issue in greater detail at the end, where we 
also show how tunable repulsive long-range interactions 
could be mediated via photons in an optical cavity. 



II. MODEL 

Let us consider an L x L square lattice with peri- 
odic boundary conditions (a 'torus'), and place a spin-| 
on each of its 2L^ edges. Starting from Kitaev's 'toric 
code^^, we consider the more general stabilizer Hamilto- 
nian 

Ho = ^'^Upp'npnpf + ^'^Vss'nsns', (1) 

pp' ss' 

where Up = (1 - n»ep ^^^.0/2, = {1 - U^es'^^,^)/'^^ 
and crx.i,a-z,i denote the usual single-spin x and z Pauli 
operators applied to spin i. The indices p and p' run over 
all 'plaquettes' (involving the four spins on the edges of 
a unit cell), whereas s and s' run over all 'stars' (involv- 
ing the four spins around a corner of a unit cell), see 
Fig. [ij The operator rip (ug) has eigenvalues 0,1 and 
counts the number of plaquette- (star-) anyons at site p 
(s). The fourfold degenerate energy levels encode two 
qubits with logical operators given by Zi = Yik&e ^z,k 
and Xi — rifeef ^x,k, * = 1, 2, where €i and are strings 
of spins topologically equivalent to single loops around 
the torus (see Fig. [ijfor an example). These operators 
commute with all rip and and obey themselves the 
usual spin commutation relations. Note that by special- 
izing to Upp' = 2J5ppi and Vss' = 2J5ss', where J > is 
the single-anyon excitation energy, Kitaev's original toric 
code model is recovered. 

Since all rip and are mutually commuting, the 
Hamiltonian Eq. ([T]) describes two independent lattice 
gases of plaquettes and stars, respectively. Without loss 
of generality, we can thus restrict our analysis to the dy- 
namics of plaquettes and their influence on one of the Zi 
operators, say Zi = Z. A corresponding logical operator 
Zee is defined by the error correction procedure (see Fig.[l] 
and Appendix |A]) . Consequently, we set Vss' = for all 
stars while assuming the plaquette interactions Uppi to 
be of the generic form 

Upp- =2J5pp, + ^(l-^ppO: (2) 
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FIG. 1: (Color online.) Quantum memory based on the toric 
code. Illustrated is an 8 x 8 lattice (periodic boundary con- 
ditions) with a total of 128 spins- 1 (grey dots) on its edges. 
The four-body plaquette and star operators are indicated in 
transparent grey and orange, respectively. A particular choice 
for all logical operators Xi, Z\, X2, and Z2 is shown, al- 
though we will focus only on the decay oi Z\ = Z (see main 
text). A number of spins is affected by cr^-errors (green and 
red dots), leading to excited plaquettes, or 'plaquette anyons' 
(solid green and red plaquettes). Measuring the plaquette 
operators yields the positions of the excited plaquettes, but 
reveals no information about how they were originally paired 
or which path (indicated by the plaquettes framed in green 
and red) they took. A minimum-weight error correction pro- 
cedure (see Appendix [A| applies (Ta;-operators to the spins 
marked by orange circles. While the green anyons are an- 
nihilated 'properly' (with a trivial loop of errors remaining 
from the top pair and no error from the bottom pair), the 
red pair is connected around a topologically non-trivial loop 
on the torus. Although the red pair is annihilated as well, an 
uncorrected cri,-error remains on the logical Z string, having 
thereby introduced a logical error in the state stored in the 
memory. 



where Tpp' denotes the shortest distance on the torus be- 
tween the centers of plaquettes p and p' , see Fig. [l] The 
strength of the repulsive plaquette interaction is given by 
the energy v4 > 0, and the interaction is long-range for 
< a < 2 (see below). The model is also equivalent to a 
long-range Ising model, see Appendix [B] 

We model the interaction of the system with a thermal 
environment by coupling each spin to a bath which can 
introduce a x-&:^ots^ \a the initial state |5'o). From a 
stan dard master equation approach in the weak coupling 
limit^^, we derive a rate equation for the probabilities 
Pm of the system to be in state Y^m) = Diem ^a;,i|*o), 
where {to} is the set of all possible patterns of (Tj;-errors. 
This rate equation reads 



[7(-w,,(m))p^.(,„) - 7(a;,(TO))p„ 



(3) 
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FIG. 2: (Color online.) Decay of the logical Z operator in a 
non-interacting toric code. The simulation data is obtained 
for grid sizes L increasing by powers of two from 16 (blue) 
to 512 (red). All curves are ensemble averages over 10* runs. 
The main plot displays (Zee), which is the average value of Z 
one would find if an error correction scheme would be applied 
at the readout time t. The inset shows the expectation value 
of the bare (uncorrected) logical Z operator. We have used 
the parameters J = 1, T — 0.3, and 7(0) — 7(2J) — 1. See 
Appendix |A] for further details on the simulation. 



where we have defined Xi{m) to be the state to with 
an additional cra;-error applied to spin i, and uJi{m) = 
Em — ^xi{m) is the energy difference between the states 
TO and Xi{m). The time evolution of the probabili- 
ties Pm determines the decay of the expectation values 

(^(cc)) EmPm(*m|^(ec)|^'m)- 

The rates 7(w) describe the transition probabilities be- 
tween states with energy difference w. A standard expres- 
sion for 7(0;) c an be obtained from a spin-boson model 
and readsPiSl 



7(0;) = 2g„ 
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(4) 



Here, /? — l/T, with T being the temperature of the bath 
(we set Boltzmann's constant to one). For simplicity, we 
assume in the following a large cut-off frequency — ?> cxd 
and choose appropriate time units for which g„ = 1. For 
n = 1, the bath is called 'Ohmic', whereas for n > 2 it is 
called 'super-Ohmic'. We find in this work that n has a 
strong influence on the decay times of the encoded states. 



III. NON-INTERACTING CASE 

We first briefly discuss the features of a non-interacting 
system, i.e., A = Q. In this case, the relevant rates en- 
tering Eq. ^ are 7(0) (rate for an anyon to hop to a 
free neighboring site), 7(— 2J) (rate to create an anyon 
pair) and 7(2 J) = 7(— 2 J)e^-^^ (rate to annihilate a pair 
of adjacent anyons, obtained from the detailed balance 
condition). The diffusive motion of the anyons causes 
the bare and error-corrected logical operators {Z) and 
{Zee), respectively, to decay with time (see Fig. [2] and 
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Appendix [C|) . We can relate this process to a diffusion 
constant D which is generally given by D = 7(0). How- 
ever, we additionally find D — 47(— 2J) due to 'indirect 
hopping' (motion of an anyon due to a neighboring pair 
creation, see Appendix [Pj) . This rate becomes relevant, 
e.g., for a super-Ohmic bath, where 7(0) = 0. 

We obtain the lifetime as follows. The error correc- 
tion fails when the fraction / of spins with cr 2: -error is 
larger than some critical value /J^. Assuming N diffus- 
ing anyons present in the system, the accumulated num- 
ber of CTa^-errors after a time t is roughly given by NDt. 
This gives / ~ NDt/2L^, from which we determine the 
lifetime r as 



2/c 



1 



max{7(0),47(-2J)}' 



(5) 



Here, we have replaced the factor N /L'^ by the Fermi 
equilibrium density (rip) = l/(e^'^ + !)• In Ref.l^ an up- 
per bound /c < 0.11 was obtained by assuming an inde- 
pendent error model. This value yields r ~ 5.8 for the 
same parameters as used in Fig. [2] Equation (5]) gives 
reasonable values of t and is in agreement with lifetimes 
obtained from alternative derivations^ . 

Note that r is independent of the system size L, con- 
sistent with previous finding^SHlS!_ This fact is also con- 
firmed by our simulations, as shown in Fig. [2] where 
{Zee) clearly approaches a step- function with increas- 
ing L. We also see that the bare expectation value 
{Z) decays even faster with larger L. Indeed, at suf- 
ficiently short times t <^ 1/ max{7(0), 47(— 2 J)}, when 
anyon pairs have not yet diffused apart from each other 
(the 'nonsplit-pair' regime, indicated by an asterisk), 
we obtain (Z) ^ (1 - l/L)^*/2 ^ e-^*/2i. By using 
iV* ~ 4L^7(— 2J)i, it follows that (Z) decays exponen- 
tially with L. 



IV. INTERACTING CASE 

We now turn to the interacting case A > 0. We first 
determine the equilibrium number of anyons N within a 
mean-field treatment (mean-field values will be indexed 
with a subscript 'mf'). This becomes accurate in the rele- 
vant limit of large L. We obtain the single-particle energy 
at plaquette p a.s Cp = SHo/Sup = J + J2p'^pUpp'''^p' ■ 
Replacing np' by the average value rimf = N^f/L'^ and 
taking the continuum limit, we find the mean- field value 
for e„ to be 



f A 

Emf = J + ?lmf / —dr ^ J + TlmfTLa, (6) 
JLxL T 

where we use the notation La — Ca/3AL'^~°' . The con- 
stant Ca is a geometrical factor of order one, given by 
the integration of 1/r" on a unit square centered at the 
origin. In particular, cq = 1. On the other hand, we 
have rimf = l/(e^'^'"' -I- 1) since the occupation numbers 



Up can only assume the values or 1. By using Eq. ^ 
to calculate nmf, we find the self-consistent equation 
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"■mf 



with the following expansion at large La 



1 



flnL. 



In In Ln 



f3J- 



(7) 



(8) 



Higher order terms in the square parenthesis are small 
if In Lq, 3> (3 J, |lnlnia|. For fixed temperature T and 
interaction strength A, these conditions are always sat- 
isfied at sufficiently large L since La oc L^^". We have 
confirmed the validity of the mean-field approximation 
by Monte Carlo simulations, see Appendix [E] 

From Eq. Q we obtain that, even though the number 
of anyons grows with the system size L, the anyon 
density rimf goes to zero for long-range repulsive interac- 
tions with < a < 2. Hence, the population of anyons is 
increasingly diluted and the system is essentially frozen 
in the ground state at large system size. This remarkable 
effect can be attributed to the divergence of the excita- 
tion energy Cmf — Tin La, which is self-consistently de- 
termined from the anyon population in the whole sample 
due to the long-range nature of the interactions. Note 
also that, despite the fact that emf is diverging, the total 
excitation energy density nnifenif/2 goes to zero for large 
L. 

Secondly, the divergence of e,„t leads to a vanishing 
anyon pair creation rate at large L, 



7(-2e„,f) ~ T 



„(21ni„)"+2 



2LI 



(9) 



This fact allows us to revise the lifetime for the non- 
interacting memory Eq. ([5|, simply by substituting J 
with the equilibrium value e^f , yielding 



max{7(0),47(-2e,„f)}' 



(10) 



From this we obtain the lifetime of an interacting memory 
in case of an Ohmic (n ~ 1) or super-Ohmic {n > 1) bath 
as 



fcLo 



TlnL, 



T^{2\nLa, 



n+3 ' 



Ohmic 



super-Ohmic 



(11) 



in the limit of large grid size [see after Eq. ([s])]. It is clear 
from these expressions that the memory lifetime is di- 
verging with L, in strong contrast to the non-interacting 
case where it was bounded by a constant. In the Ohmic 
case, this divergence of r is entirely due to the vanishing 
density, since 7(0) — 2T is non-zero. In the super-Ohmic 
case, however, an additional divergence due to the vanish- 
ing of 7(— 2enit) is obtained, see Eq. Since the energy 
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gap grows logarithmically with L, t grows polynomially, 
but with a rather favorable power. For instance, constant 
interaction (a = 0, see also below) leads to t oc L? /\nL 
in the Ohmic case and to r oc L^/ln^L in the super- 
Ohmic [n = 2) case. 

Eq. (10) is valid in the nican-ficld limit and does not 
include fluctuations which give rise to additional errors 
due to the long-range interactions. However, the mem- 
ory remains self-correcting both for an Ohmic and for a 
super-Ohmic bath. Indeed, for an Ohmic bath, we can 
neglect the effect of the repulsive force if the change of 
energy w in a diffusive step is smaller in magnitude than 
T (see Eq. Q), so that we can approximate 7(0;) ~ 7(0). 
In particular, for a single pair of anyons at distance r, we 
have |a;| < aA/r"''"^, which defines a critical radius 
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(12) 



beyond which the fluctuations become negligible. For 
a = one has Tc = 0. For a > 0, since the average 
distance ~ 1/ ^nmf between anyons grows with L while 
is independent of L, the fluctuations also become neg- 
ligible. The validity of Eq. (10) for the Ohmic case is 



confirmed by numerical simulations both for a = (see 
Fig. |3]) and a > (see Fig. [I]). 

Concerning the super-Ohmic case, we expect devia- 
tions from Eq. (10) if the fluctuations of a; ~ are more 
effective for the anyon motion than the indirect diffusion 
mechanism, proportional to the rate ([9|. However, due 
to the decreasing interaction strength, such fluctuations 
in Lj become small at large L and still result in a van- 
ishing diffusion coefficient. Furthermore, at a = direct 
hopping is impossible and Eq. (10) is valid (see Fig. [s]). 
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FIG. 3: (Color online.) Thermal stability of the interacting 
memory. The data in the top (bottom) panel was obtained 
for an Ohmic (super-Ohmic, n = 2) bath. Plotted as a func- 
tion of L are the numerically simulated times at which the 
expectation values of the bare (squares) and error-corrected 
(diamonds) logical Z operator have decayed from 1 to 0.9. 
The dotted lines serve as a guide to the eye. The red dashed- 
dotted curves are calculated from Eq. (10 1 with /c = 0.11, 
where we have used the self-consistent values of n^f and 
from Eqs. (|6| an d |7[ ). Similarly, the green dashed lines are 
also due to Eq. ( |10[ |, but here /c was fit to the numerical 
data of the 90% threshold times, yielding /c = 0.022 for an 
Ohmic, and fc = 0.007 for a super-Ohmic bath. The inset 
shows the decay of (Zee) with time, and the 90% threshold 
is illustrated by the dotted line. It is seen that choosing this 
particular value has no substantial influence on the scaling 
behavior with L. Parameters used in these simulations were 
J = IA= 0.1, and r = 0.3. 



V. NUMERICAL SIMULATIONS 

We turn now to the numerical simulations of our 
model, Eq. (IT]), and focus first on constant long-range 
interactions (a = 0). In this case, the total energy 
E]s — N J + ^N{N — 1) depends only on the number 
of anyons N , but not on their position. This simplifies 
the numerical treatment considerably. Our results are 
displayed in Fig. |3] The numerical data show a clear in- 
crease of the memory lifetime r with L. Note that this 
holds already for the hare Z . Like in the non-interacting 
case (see Fig.[2|, the beneficiary effect of the error correc- 
tion at read-out is to prolong the lifetime by maintaining 
{Zee) close to one (see inset of Fig. [s]). 

Our analytical results describe the numerical data re- 
markably well. By fitting fc in Eq. ([lO| to the simu- obtain T7(-2J) ~ 5 x iC^, i.e., already for a moder- 
ate system size the lifetime t of the memory is about 
a million times longer than the single-spin lifetime. For 
quantum dots, the latter is typically in the range of mil- 
liseconds to seconds at about 100 mKl^^E^. 

For non-constant long-range interaction (0 < a < 2), 
simulating the time dynamics of the memory is numeri- 
cally more costly due to an 0{LP') overhead coming from 
recalculating all spin flip rates. Nevertheless, we were 



The lifetime t can be compared to the physical time 
scales of single spin flips, 1/7(0) and l/7(— 2J). For in- 
stance, for the L — 256 super-Ohmic case in Fig. 3 we 



lation data, excellent agreement is found for an Ohmic 
bath (top panel of Fig. |3|, while for a super-Ohmic 
bath (lower panel), analytics and numerics agree well for 
L > 64. Further, the fit yields values for fc of about 
0.01 — 0.02, which is reasonable in comparison to the 
upper bound fc = 0.11 found for a model of uncorre- 
lated errors (dashed-dotted lines in Fig. [sjP'. See also 
Appendix [F] for an extended discussion. 
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FIG. 4: Thermal stability of the interacting memory 
with a 7^ and an Ohmic bath. Data points refer to 
the numerically calculated times at which the error-corrected 
logical Z operator has decayed from 1 to 0.9 in the cases 
a = (diamonds), a = 0.5 (triangles), and a — 1 (squares). 
Note that we have replotted the data from a = merely 
for comparison. The dashed lines are due to Eq. 1 10 1 (as in 
Fig. [3|, with a fit of fc yielding fc = 0.027 for a = 0.5 and 
fc ~ 0.032 for a — 1. Inset: Decay of (Z^c) as a function of 
time for different grid sizes, L = 8, 16, . . . , 256, and a = 1. 
Parameters used in the simulations were J = 1, A = 0.1, and 
r = 0.3. 



able to study the (more tractable) case of an Ohmic 
bath. The results are presented in Fig. [4] for a = 0.5 
and a — I. Clearly, the memory lifetime is still increas- 
ing with L, proving the memory to be self-correcting also 
for a 0. Furthermore, the data is in perfect agreement 
with the analytically calculated lifetime Eq. (10). The 



super-Ohmic case for a > is more difhcult due to the 
increased memory lifetime, and will be examined else- 
where. 

We consider now in greater detail the super-Ohmic case 
at a = 0, which has the most favorable scaling. On a 
sufficiently short time scale (in the nonspUt-pair regime), 
the rate equation 



dN, 



mi 



dt 



-24t)"^mf7(2e,;f), 



(13) 



describes the initial time-evolution of the system well, 
since in this non-diffusive regime only pair creatiorfi^l and 
annihilation takes place. We confirm this by comparing 
a numerical integration of Eq. ( 13 ) with a direct simula- 



tion, presented in Fig. [5] After a rapid initial 'build-up' 
phase, N*^f saturates to a value determined by the self- 
consistent condition N*^f = 4L^e~^('^+'^^»f)^, obtained 
by setting dN^Jdt — in Eq. (13). In this state, 
the excitation energy is diverging with L, since we have 
e*^f ~ AN^f ~ AN,^f/2 cx \nL. This effectively sup- 
presses the indirect diffusion of anyons. Therefore, the 
system remains in a quasi-stationary state which evolves 
to the final anyon density on a time scale also diverg- 
ing with L. In this regime of nonsplit pairs, one has 
(Z) ~ e~^">f/^^. This leads to the quasi-stationary value 



Fig.[5| 
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which approaches one for large L (see 



VI. DISCUSSION AND CONCLUSIONS 

We now discuss various physical aspects arising in the 
context of our model. We first remark that super-Ohmic 
baths, which provide the best scaling of the lifetime with 
L, are not uncommon and emerge, e.g., for quantum dot 
spins in contact with phonona^. As for the periodic 
boundary conditions, these are no t an essential ingredi- 
ent to a topological stabilizer cod d^^ l ^H Concerning the 
many-body nature of the interactions involved, general n- 
body couplings can in pri nciple be engineered from short- 
range 2-body interactiona^SnSni Por e xamp le, toric codes 
with interacting anyons are derived iiPI^H. A systematic 
procedure to construct such effective low-energy Hamil- 
tonians can be rigor ously founded on the Schrieffer- Wolff 
transformatior(2212il. A promising candidate system to re- 
alize topological models are ultracold atoms or molecules 
in optical lattices^^. 

In the same way, physical long-range interactions of 
the type considered in this work could also be gener- 
ated perturbatively. A well-known example is the RKKY 
interactiorP^, e.g., for a 2D Kondo-lattice of nuclear 
spinJ^. Alternatively, constant interactions (a = 0) 
can be realized for qubits coupled to photon modes in 
QED-cavitie j^ " i The interaction range is determined 
by the wavelength of the photon and can reach macro- 
scopic distances, in particular in superconducting cavity 
stripline422H25|. As an illustrative example let us consider 
a two-photon coupling (see also Appendixlcl) of the form 



2 

^^int ^^i^iajai + ^gpnp{ala2 + ai4). 



(14) 



Here, U!i are the photon frequencies, and gp is the cou- 
pling strength of plaquette p to the modes. Note that Or, 
is spatially constant over the photon wavele ngth Xp^. 
Performing a Schrieffer- Wolff transformatiorP^ESl (gee 
also Appendix [C]) yields an effective Hamiltonian con- 
taining the desired constant anyon interaction. 



Hctf — 



a\ai 



UJl — UJ2 



2 

E 



ujia\ai 



(15) 



The strength and sign of the interaction is tunable via 
the difference in frequencies and occupation numbers of 
the modes, and can consequently be made repulsive in 
a steady-state regime. For typical stripline cavities with 
\i ~ cm and for typical lattice constants of 100 nm (e.g. 
quantum dots), we see that the anyon interaction stays 
constant over system sizes L as large as 10^. With the 
super-Ohmic scaling (r ~ L^), and a single-spin lifetime 
l/7(— 2J) ~ 1/is — 1 ji8|36 | ^]^Q memory lifetime r can be 
made sufficiently long for any practical purposes. 
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FIG. 5: (Color online.) Short-time dynamics of the interact- 
ing memory in a super-Ohmic bath. In this case, the mem- 
ory is in the nonspht-pair regime. The curves refer to differ- 
ent values of L increasing from L = QA (blue) to L = 2048 
(red) in powers of two. Upper panel: The time dependence 
of the anyon number A'^ obtained from the simulations (sohd 
lines) is compared to the solutions of Eq. (131 (dashed hues). 
The crosses are the exact values A'^* obtained from a par- 
tition function of pairs (see Appendix [e|). Good agreement 
with A" is also obtained for the lower curves at longer times 
(not shown). Lower panel: The expectation value of the bare 
Z obtained from the simulations (solid lines) is compared to 
g-jv /2L ((;jagijg(j lines), where N*(t) is obtained from the up- 
per plot. Parameters used are J = 1, A = 0.1, and T = 0.3. 



In conclusion, we have demonstrated the existence of 
2D stabihzer quantum memories at finite temperatures. 
The stability of the memory results from the repulsive 
long-range nature of the interaction we have introduced. 
We expect that similar systems in the presence of such 
interactions also prove useful as self-correcting quantum 
memories. 
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Appendix A: Numerical simulation 

1. Time evolution 

In order to achieve a time evolution in accordance 
with Eq ([3]), each iteration of a simulation consists of 
the following steps, (i) We record the relevant parame- 
ters of the system, (ii) We calculate the total spin flip 
rate R = J2i^i^s ~ ^xi(s))j where s is the current state 
of the system, (iii) We draw the time At it takes for 
the next spin to flip from an exponential distribution, 
At ^ Exp(l/i?), and then add this to the current total 
time, (iv) We calculate all individual spin flip proba- 
bilities Pi = 7(es — ej;.(^s))/R and flip a spin at random 
accordingly. After some initially specified time has been 
reached, we stop and have obtained a single 'run'. The 
final data presented in this work is then generated by 
averaging over many (typically several thousand) runs. 

2. Error correction 

We have seen that, for an interacting system, already 
the bare logical Z operator becomes stable with increas- 
ing L. Nevertheless, it is useful to apply an error cor- 
rection scheme once the memory is being read out. By 
(Zec)(i), we denote in this work the average value of Z 
we would have obtained if we had performed error cor- 
rection at time t. The goal here is to properly annihi- 
late corresponding anyons (by applying cr ^^-operations), 
thereby reverting the undesired operations performed by 
anyon paths crossing the logical operator strings. How- 
ever, since only the positions of the anyons are known, 
this correspondence has to be guessed. We do this by 
choosing the pairing with the minimal sum of connec- 
tion path lengths using Blossom which is the lat- 
est improvement on Edmonds' minimal-weight perfect 
matching algorithm^. If many anyons are present, us- 
ing the complete graph as the input to this algorithm 
is numerically infeasible. In excellent approximation, 
we therefore replace the complete graph by a Delaunay 
triangulation^Si. 



Appendix B: Mapping from lattice gas to Ising 
model 

Note that Hq in Eq. ([T]) has the general form of two 
independent lattice gases, which is also equivalent to two 
Ising spin lattices. We explicitly perform the transfor- 
mation in the plaquette sector by identifying the Ising 
variables Sp = I — 2np, yielding 



Ho 



p' 



p,p' 



. . . . . ^^^^ 

where Upp' is given in Eq. ([2]) and the primes in the sum- 
mations indicate p' 7^ p. We have used Upp = 2J and 
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Up' p. The noninteracting Kitaev model corre- 



sponds to noninteracting spins in an external magnetic 
field. The ground state corresponds to the fully polar- 
ized state Sp — 1 for all p, where no anyon is present. 
However, for T > a finite density of anyons emerges 
and is sufficient to destroy the information stored in the 
memory. 

If a short-range ferromagnetic interaction is intro- 
duced, ordering of the system is spontaneously favored 
below some critical temperature. A higher magnetization 
corresponds to a lower population of anyons and improves 
the lifetime. However, short range interactions do not 
improve the scaling of the lifetime with the system size, 
since a residual density of anyons is left at any finite tem- 
perature. As in the noninteracting finite density 
of excited plaquettes efhciently destroys the stored quan- 
tum information, in agreement with the general analysis 
Instead, repulsive long-range interactions lead to 
a fully polarized system at a given temperature for suffi- 
ciently large system size L. 



Appendix C: Lifetime in the presence of a single pair 

The decay of the bare and logical Z operators is most 
simply illustrated by assuming only a single anyon pair in 
the memory. We set 7(2 J) = 0, so that pair creation and 
annihilation are not allowed. If no anyons were present, 
the initial values (Z) = (Z^c) — 1 would be stable. We 
apply one tTa;-operation at a randomly chosen site and 
thereby create two neighboring anyons at t = 0. This 
causes a partial decay of the bare logical operator already 
at t = 0, since we might have chosen to flip a spin on the 
logical Z operator, yielding (Z) = 1 — This has been 
used in the main text in the discussion of the nonsplit- 
pair regime. 

We now study the decay for t > in the continuum 
limit and therefore neglect the 1/L correction at t — 0. 
We consider a single pair of diffusing anyons with coordi- 
nates {xi,yi) and (a;2, j/2) created at the origin. We then 
assume that the probability to find an anyon at position r 
is described by the probability density 



p{r) 



1 

47r7(0)t^ 



(CI) 



We represent the torus as an infinite plane with the 
points {x,y) and {x + mL^y + nL) being equivalent 
{m,n g Z). The logical Z operator is then represented 
by parallel lines at yz — L/2 + nL. The two anyons 
diffuse along y with probability density p{y.i — j/o) = 
e-fa.-yo)V47(o)t/^47r-y(o)t, where i = 1,2 and the ini- 
tial (random) coordinate satisfies —L/2 < yo < L/2. 
The average of the logical operator at time t is 



(Z) 



f-L/2 



-L/2 



dyo 
L 



dyidy2p{vi - yo)p{y2 



yo)z{yi,y2), 
(C2) 



A 0.5 




0.05 0.1 

j{0)t/L^ 



0.2 



FIG. 6: (Color online.) Decay of the bare and corrected ex- 
pectation value of Z due to a single pair of anyons in the mem- 
ory. The dots show numerical data (averaged over 10^ sam- 
ples) while the two curves are the continuum limit expressions 



Eq. (|C3| and for (Z^c) (solid) and (Z) (dashed). The nu- 
merical data has been obtained for L = 32 with 7(0) = 1, 5, 10 
and L = 64, 128 with 7(0) = 1. All points collapse onto each 
other when plotted as a function of 7(0)t/Z/^. 



where z{yi, j/2) gives the sign of Z if the two anyons have 
diffused to the coordinates yi and 2/2- Since Z changes 
sign each time an anyon crosses the lines at yz, we have 
z{yi,y2) = z(yi)z{y2) where z{y) = 1 if -L/2 + 2nL < 
y < L/2 + 2nL and —1 otherwise (71 G Z). Therefore we 
can write 

.1/2 

{Z) = / dzofizn)^, 
J-1/2 



where we have made the change of variables yo 
such that 



(C3) 
Lzq, 



fizo) 



+ 00 



E (-1)" 



erf 



erf 



2zo -I- 2n + 1 
A^j{0)t/L\ 



2zn 



2n-l 



4^7(0)^/^2 



(C4) 



We now consider the average of the error-corrected 
logical operator Zee- In this case, only the distance 
2/12 — Vi ~ y2 between the two anyons is important since 
the value of Z^c is 1 if ~L/2 + 2nL < yu < L/2 + 2nL, 
and is —1 otherwise. The probability distribution for yi2 
is / dy2piyi2 - y2)p{y2) = £-^'12/8^(0)7^8^7(0)^, which 
gives 



+00 



(Ze. 



2n-t- 1 



2^27(0)7^2] 



(C5) 



Both functions (C3| and (C5| are plotted in Fig. |6] and 
show perfect agreement with the numerical simulation. 
An important feature of the above analytical expressions 
is that the time dependence only enters through the com- 
bination j{0)t/L'^, which makes it possible to scale curves 
from different system sizes and diffusion constants onto 
each other. 
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Appendix D: Indirect diffusion of anyons 

To determine the rate D associated with indirect dif- 
fusion, we consider an isolated anyon in the lattice and 
its probability pij to be at site In the Ohmic case, 

we have 7(0) 0, and direct hopping to neighboring 
sites is thus allowed. In the continuum limit, a stan- 
dard diffusion equation '^^^^ = DV^p(r) with D — 7(0) 
is obtained. In the non-interacting super-Ohmic case, 
however, 7(0) = 0. We assume 2/3J ^ 1, such that. 



since 7 (2 J) 



7(— 2J), the recombination of a pair 



of anyons is essentially instantaneous. Hopping from the 
site to, e.g., + 2) is possible by creation of an 
anyon pair occupying sites + and {i,j+2), an event 
occurring with rate 7(— 2J). Since the intermediate state 
can decay back to the initial state, the actual rate for the 
indirect hopping process is 7(— 2J)/2. Similar consider- 
ations hold for all other sites. Accounting for all these, 
we write 



dpi 



7(-2J) 



dt 



(-12p, 



-Pi+2,. 



+2pi+ij+i + 2p, 



Pi,j+2 ^ 

1 + 2p,_ 



P. 



which in the continuum limit yields D = 4q{—2J). We 
can expect that the properties of the memory improve 
by lowering the value of 7(0), but only as long as 7(0) > 
47(— 2J). In the interacting case, J is replaced by an 
appropriate excitation energy (e.g., e^i at equilibrium). 



Appendix E: Equilibrium density of interacting 
anyons 

The equilibrium number of excited plaquettes can 
be approximated with arbitrary accuracy by using the 
Metropolis algorithm^ to sample the probability dis- 
tribution cx 6-/^/2 Ep.p't/pp'^prv^ see Eq. This can 
be used to study the accuracy of the mean-field value 
-^mf = n^iL? [see Eq. ([t])], in particular for values a ^ 0. 

Due to the long-range nature of the interaction, we find 
that iVmf compares very well to the equilibrium value of 
N at generic values of the temperature and interaction 
exponent a. This is illustrated in Fig. [7j which shows a 
satisfactory agreement already at moderate values of L. 

We also note that, for the case of constant interaction 
(a = 0), the average number can be calculated directly 
from the grand-canonical partition function 



E 

2fc<L2 



L2 

2k 



(El) 



since the energy of a given anyon configuration does not 
depend on the positions of the anyons, but only on their 
total number N — X]p"p- presence of an ap- 

preciable anyon interaction and at low temperature, the 
number of excited plaquettes is much smaller than L^. 
Therefore, one can restrict the sum (El ) to the first rel- 
evant terms. 
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FIG. 7: (Color online.) Comparison of the equilibrium value 
of A'^ obtained numerically (crosses) with A'^mf (curves) for 
different grid sizes. We have used the interaction exponents 
a = (solid line), a = 0.5 (dashed line), and a — 1.0 (dotted 
line). Other parameters are J — 1 and T = 0.5. 



In a similar way, the exact quasi-stationary number of 
paired anyons N* , displayed in Fig. [s] (crosses), can be 
calculated from a partition function 



E 

fc<2L2 



2L^ 



'I3E2 



(E2) 



Here we have assumed that k sufficiently diluted errors 
are present in the memory such that 2k anyons are cre- 
ated in the nonsplit-pair regime. The average number of 
anyons N* calculated from (E2) is in very good agree- 
ment with the simulations, see Fig. [5] 
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FIG. 8: (Color online.) Average of the corrected operator 
^ec for a model with independent (Ji,-errors occurring with 
probability / at each spin. The solid curves refer to our nu- 
merical simulations with lattice sizes L = 40, 100, 200. The 
error correction fails at a value fc — 0.1, which is shghtly 
smaller than the value 0.11 from Ref.^ (dotted line). In the 
inset, we plot the value of fc from simulations of the non- 
interacting toric code in contact with a bath at temperature 
T (7(0) = 7(2 J) = 1). The fraction fc is extracted at the time 
r when {Zee) decays to zero in the limit of large L (see Fig.[2|. 
This value is always smaller than / = 0.11 and depends on T. 



Appendix F: Critical fraction of errors 

We discuss here in greater detail the role of fc in the 
memory lifetime Eq. ([5| . An analogous result can be ob- 
tained based on the following different reasoninjl^l, xhe 
distance between the two anyons of a pair after a time r 
is of order A£ = \/ Dt and this is required to be much 
smaller than the average anyon separation ~ L'^/N. 
This gives 



T < 



max{7(0),47(-2J)}' 



(Fl) 



Interestingly, this upper bound coincides with the right- 
hand side of Eq. ^ if the probabihty for each spin to be 
flipped is fc = \- Hence, this should be an upper bound 
on the values of fc we determine numerically. In the 
following, we demonstrate that this is indeed the case, 
and that the values we obtain are meaningful. 

We first determine fc for an independent error model. 
By mapping the toric code to a random-bond Ising 
model, a phase transition at fc — 0.11 was obtained 
in Ref.'^. Numerically, we find in this case the value 
fc ~ 0.1, see Fig. |8] This shows that our minimum- 
weight error correction scheme works close to optimal. 

In the simulations of the non-interacting model, we 
observe a sharp transition in time similar to Fig. [s] (see 
Fig. [2]) but, not surprisingly, we find at the transition 
point a relative fraction of errors different from 0.1. We 
attribute this fact to the statistical properties of the 
distribution of errors generated by the time evolution. 
Clearly, the errors created by the anyons in their diffu- 
sive motion have strong spatial correlations, rather than 



100 




10 r 



T 



FIG. 9: The values of r extracted at the sharp transitions 
of the {Zac) decay (circles) are compared to Eq. ([5| (dashed 
curve). Very good agreement is obtained for fc — 0. iP. 



being independent and uniformly distributed across the 
sample. We find that such correlations yield values of fc 
strictly smaller than 0.1 but still of the order of a few 
percent, see Fig. |8] 

We mention here a second uncertainty in Eq. ([5|, 
namely the estimate NDt of the total number of errors. 
An isolated anyon can have either one or three tr^-errors 
at its plaquette spins. In the first case it contributes to 
the error rate with 2D = 3D — D since for three hopping 
processes one error is created and for the fourth a pre- 
existing error is removed. If three (Ta;-errors are present, 
an opposite rate —2D is obtained but three-error plaque- 
ttes are rare in the initial time evolution (they only arise 
from the crossing of two anyon paths). This gives a to- 
tal error rate from diffusion of order D. Furthermore, a 
typical number of anyons is determined by the equi- 
librium density, but the actual value is always smaller 
(since there are initially no anyons present in the mem- 
ory). Therefore, NDt 
number of errors by an unknown prefactor 



j^^j differs form the actual 



From the previous discussion it is clear that fc in 
Eq. ([5]) is generally different from the value fc = 0.11 
known in the literatur^. This justifies our practice of 
using fc as a fitting parameter to reproduce the func- 
tional dependence of the lifetime, e.g., as a function of 
L or T. Remarkably, the asymptotic dependence on L 
is well described by Eq. ([5]), both in the interacting and 
non- interacting case (see Figs. [2] and [s]). An example of 
the temperature dependence of t in the non-interacting 
case is shown in Fig. |9] and is also described by Eq. ^ 
very well. Note that fitting the data always yield values 
of fc smaller than fc — 0.11, and of the order of a few 
percent. This is consistent with their interpretation as a 
critical fraction of errors. 
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Appendix G: Effective Hamiltonian via 
Schrieffer- Wolff transformation 



In order to find an effective Hamiltonian for (14), 



we write H = Hq + V , where Hq = Y^i=i ^iO^iO-i and 
V = '^pgpnp{a\a2 + a,ia\), and treat F as a small per- 
turbation. The general expression for the Schrieffer- Wolff 
transformation of H up to second order in V reads 



Ho 



lim 



dte-''[V,V{t)]+0{V^), (Gl) 



where V{t) = exp{iHQt)V exp{—iHQt), which yields in 
our case 

V{t) =Y,9pnp (e*("i-"^)*ala2 + e-'("i-'^^)*4ai) . 



With this, the commutator in Eq. (Gl ) evaluates to 



(G2) 



[V, V(t)] ~ 2iC^ gpUp)^ {ala2 ~ a\ai) sin(a;i — uj2)t. 

(G3) 



Inserting this into Eq. (Gl) and performing the integral 
yields Eq. (flSl). 



A two-photon coupling of the type a\a2 can emerge 
from a quadratic coupling of spin to electric (or magnetic) 
cavity fields such as E^Ey. Now, the toric code model 
can be obtained from Kitaev's honeycomb model by per- 
turbation theorjBi. Indeed, in leading order one obtains 
J oc J'^Jy/Jl, where Jk are the exchange couplings in the 
honeycomb modeP'^. In multiferroic materials electric 
fields can couple to the spin (-texture) via a modification 
of the exchange interaction like Jk ^ Jk + IkE]^ (with 
lx,y,z some coupling constants). Thus, if e.g. one Jx and 
one Jy occurring in J oc J^Jyl Jz get modified in this way 
(by locally modifying the corresponding links) , we end up 
with a coupling of the desired form. By going to reso- 
nance uji « CJ2, the quadratic term E^Ey can be made 
dominant over the linear ones (which are non- resonant). 
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